function friction

%  solve IVP using MATLAB routines
%         y' = f(t,y)  with   y(0) = y0        '
%  where y = (y1, y2, y3 , ..., yn) is an n-vector

clear *
clf

% time interval
tmax=1000;

%  initial values
y10=1; y20=1; y30=0.1;

%  calculate solution using a MATLAB routine
[t,y] = ode45(@rhs,[0 tmax],[y10 y20 y30]);
%[t,y] = ode23s(@rhs,[0 tmax],[y10 y20 y30]);

subplot(2,1,1);
hold on
plot(t,y(:,1),'r')
box on
% commands to label each axes
xlabel('t-axis','FontSize',14,'FontWeight','bold')
ylabel('y-axis','FontSize',14,'FontWeight','bold')
% command to put legend into plot
%legend(' y_1',' y_2',1);
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 
% Set legend font to 14/bold                            		
%set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold'); 
hold off

subplot(2,1,2);
hold on
plot(t,y(:,3),'--b')
box on
% commands to label each axes
xlabel('t-axis','FontSize',14,'FontWeight','bold')
ylabel('\phi-axis','FontSize',14,'FontWeight','bold')
% command to put legend into plot
%legend(' y_1',' y_2',1);
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 
% Set legend font to 14/bold                            		
%set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold'); 
hold off

%  define RHS of equations
function dy=rhs(t,y)
dy=zeros(3,1);
alpha=4; delta=1; kappa=1; gamma=0.5; beta=0.25*(1+sqrt(5))*1.5;
dy(1) = y(2);
dy(2) = t - y(1) - kappa - alpha*y(3) - delta*y(2);
dy(3) = y(3)*(1 - y(3) - gamma*y(2))/beta;














